Introduction to Bayesian Disease Measurement for Health Scientists - Day 1 | 2

CA18208 HARMONY Larissa Training School 2024 - https://harmony-net.eu/

Eleftherios Meletis
Polychronis Kostoulas
Evangelia Anyfanti

2024-04-22 & 2024-04-23

Bayesian Statistics

In this session we’ll see how we can estimate a probability of interest but in a Bayesian framework, i.e. using Bayes theorem.

Bayes’ theorem - P(A|B) = P(B|A)*P(A)/P(B)

Components

P(A|B): Prob of event A occurring given that B is true - Posterior probability
P(B|A): Prob of event B occurring given that A is true - Likelihood ~ function of A
P(A): Prob of event A occurring - Prior probability
P(B): Prob of event B occurring - Marginal probability ~ sum over all possible values of A

What we usually see/use

theta: parameter of interest | y: observed data}
P(theta|y) = P(y|theta) * P(theta)/P(y)

Where:

P(theta): Prior probability of parameter(s) of interest;
P(y|theta): Likelihood of the data given the parameters value(s)
P(theta|y): Posterior probability of parameter(s) of interest given the data and the prior

Bayesian Inference - Summary & Example

To estimate the posterior distribution P(theta|y) we need to:

Specify the Prior distribution: P(theta)
Define the Likelihood of the data: P(y|theta)

Example 1

Hemophilia is a disease that exhibits X-chromosome-linked recessive inheritance. Consider a woman who has an affected brother, which implies that her mother must be a carrier of the hemophilia gene. We are also told that her father is not affected. Suppose she has two sons, neither of whom is affected.

The unknown quantity of interest is the state of the woman.

Exercise 1: Bayesian apparent prevalence (ap) estimation

y out of n individuals test positive. Estimate the apparent prevalence in a Bayesian framework.

Parameter of interest: ap - [0,1]
Data: n tested, y positive
Prior distribution for ap: ap ~ Beta(a,b)
Likelihood: y ~ Binomial(n,ap)

Let’s write our first JAGS model

ap_model <- 
  'model {
  
  # Define likelihood distribution of the data
  # JAGS Binomial distribution Arguments: p, n 
  
  y1 ~ dbin(ap1,n1)
  y2 ~ dbin(ap2,n2)
  
  # Specify prior distribution for par of interest 
  # Uniform (non-informative) prior distribution 
  ap1 ~ dbeta(1, 100000)
  ap2 ~ dbeta(1, 100000)


  #data# n1, y1, n2, y2
  #monitor# ap1, ap2
  #inits# ap1
  }
  '

Let’s run our first JAGS model

# Call JAGS
library(runjags)

# Provide Data 
n1 = 4000
y1 = 1200

n2 = 40
y2 = 12

# Initial values for par of interest

ap1 <- list(chain1=0.05, chain2=0.95)


# Run the model
results <- run.jags(ap_model, n.chains=1, 
                    burnin=5000, sample=10000)
## Warning: Convergence cannot be assessed with only 1 chain
# View results

# Plot results
plot(results)

# Print results
summary(results)
##          Lower95       Median      Upper95         Mean           SD Mode
## ap1 1.091972e-02 0.0115451934 0.0122377064 0.0115501328 3.355818e-04   NA
## ap2 6.553641e-05 0.0001266365 0.0002000996 0.0001296519 3.512836e-05   NA
##            MCerr MC%ofSD SSeff        AC.10 psrf
## ap1 4.270690e-06     1.3  6174  0.001341437   NA
## ap2 4.682192e-07     1.3  5629 -0.002783184   NA
# Convert the results to 'mcmc.list' object
mcmc_ap <- as.mcmc.list(results)

# Convert 'mcmc.list' to data frame for ggplot
df_ap <- as.data.frame(as.mcmc(mcmc_ap))

library(ggplot2)
library(tidyr)  # For pivot_longer

# Assuming df_ap contains the columns for ap1 and ap2
# If your data frame isn't set up this way, adjust the data extraction part accordingly.

# Convert the dataframe to long format
df_long <- pivot_longer(df_ap, cols = c(ap1, ap2), names_to = "parameter", values_to = "value")

library(ggplot2)

# Assuming df_long is already created and contains the parameters 'ap1' and 'ap2' in long format

# Generate a sample from a beta distribution
beta_data <- data.frame(value = rbeta(10000, 1, 10000))
beta_data$parameter <- "Beta Distribution"

# Combine this with your existing data
df_long <- rbind(df_long, beta_data)

# Plotting
p <- ggplot(df_long, aes(x = value, fill = parameter)) +
  geom_density(alpha = 0.6, adjust = 1) +  # Adjust can be tweaked for smoother curves
  scale_fill_manual(values = c(ap1 = "blue", ap2 = "red", `Beta Distribution` = "green")) +
  labs(title = "Overlay of Density Plots for Apparent Prevalences ap1, ap2, and Beta Distribution",
       x = "Prevalence Estimate",
       y = "Density") +
  theme_minimal()

# Print the plot
print(p)

Exercise 1 - Practical

What if you tested 4000 individuals (n=4000) and 1200 tested positive (y=1200)?
What if you have prior information on ap?

Exercise 2: Bayesian true prevalence (tp) estimation

Assuming the absence of a perfect test we do not know how many individuals are truly positive/negative.

Instead we know that n individuals are tested with an imperfect test and y have a positive result.

Apparent/True prevalence: ap/tp
Sensitivity: Se
Specificity: Sp
ap = P(T+) = P(T + & D+) + P(T+ & D-) = P(D+) * P(T+|D+) + P(D-) * P(T+|D-) =>
ap = tp * Se + (1 - tp) * (1 - Sp)

Create a JAGS model for true prevalence estimation

Parameters of interest - tp, Se, Sp

Prior distributions

tp ~ dbeta(1,1) # Uniform (non-informative) prior distribution
Se ~ dbeta(25.4, 3.4) # 0.85 (0.7 - 0.95)
Sp ~ dbeta(95, 5) # 0.77 (0.49 - 0.96)

Data: n tested, y positive

Likelihood: y ~ Binomial(n,ap),
ap = tp * Se + (1 - tp) * (1 - Sp)

Write JAGS model

tp_model <- 
  'model {
  y ~ dbin(ap,n)
  ap <- tp*Se + (1-tp)*(1-Sp)
  
  # Uniform (non-informative) prior distribution 
  tp ~ dbeta(1,1)
  # Informative priors for Se and Sp
  Se ~ dbeta(25.4, 3.4)
  Sp ~ dbeta(95, 5)
  
  #data# n, y
  #monitor# tp, Se, Sp
  #inits# tp, Se, Sp
  }
  '

Let’s run our JAGS model

# Call JAGS
library(runjags)
# Provide Data 
n = 4072
y = 1210
# Initial values for pars of interest
tp <- list(chain1=0.05, chain2=0.95)
Se <- list(chain1=0.05, chain2=0.95)
Sp <- list(chain1=0.05, chain2=0.95)

# Run the model
results <- run.jags(tp_model, n.chains=2, 
                    burnin=5000, sample=10000)
## View results
# Plot results
plot(results) # Click backwards to view all plots

# Print results
summary(results)
##      Lower95    Median   Upper95      Mean         SD Mode        MCerr MC%ofSD
## tp 0.2403075 0.2978380 0.3684788 0.3001641 0.03217121   NA 0.0012618569     3.9
## Se 0.7566369 0.8866616 0.9812408 0.8770171 0.06178538   NA 0.0022247881     3.6
## Sp 0.9076471 0.9524757 0.9866954 0.9497712 0.02142762   NA 0.0006325174     3.0
##    SSeff     AC.10     psrf
## tp   650 0.5536472 1.014044
## Se   771 0.4485606 1.002960
## Sp  1148 0.3088396 1.003026

Exercise 2 - Practical

What if you had no prior information at all on Se,Sp?
What if the sample size is n=40, y=12?

Ok for today? - Any questions?

Sensitivity - Specificity estimation without a gold standard

Hui-Walter paradigm/model (1980)
A particular model formulation that was originally designed for evaluating diagnostic tests in the absence of a gold standard
Not originally/necessarily Bayesian - implemented using Maximum Likelihood
Evaluating an imperfect test against another imperfect test
If we don’t know the true disease status, how can we estimate sensitivity or specificity for either test?

Hui-Walter explained

https://www.youtube.com/watch?v=z6devQmW2xE&ab_channel=PolychronisKostoulas#

We will use the data/observations from the manuscript published back in 1980.

Hui-Walter (1980) dataset Table 1

pop_1 = matrix(nrow=3,ncol=3)
rownames(pop_1) = c("Mantoux_Test_Pos", "Mantoux_Test_Neg", "Total")
colnames(pop_1) = c("Tine_Test_Pos", "Tine_Test_Neg", "Total")

pop_1[1,1] = 14
pop_1[1,2] = 4
pop_1[2,1] = 9
pop_1[2,2] = 528
#Total rows and columns
pop_1[1,3] = pop_1[1,1] + pop_1[1,2]
pop_1[2,3] = pop_1[2,1] + pop_1[2,2]
pop_1[3,1] = pop_1[1,1] + pop_1[2,1]
pop_1[3,2] = pop_1[1,2] + pop_1[2,2]
N_1 = sum(pop_1[1,1] + pop_1[1,2] + pop_1[2,1] + pop_1[2,2])
pop_1[3,3] = N_1
pop_1
##                  Tine_Test_Pos Tine_Test_Neg Total
## Mantoux_Test_Pos            14             4    18
## Mantoux_Test_Neg             9           528   537
## Total                       23           532   555

## Now let's do pop_2
pop_2 = matrix(nrow=3,ncol=3)
rownames(pop_2) = c("Mantoux_Test_Pos", "Mantoux_Test_Neg", "Total")
colnames(pop_2) = c("Tine_Test_Pos", "Tine_Test_Neg", "Total")

pop_2[1,1] = 887
pop_2[1,2] = 31
pop_2[2,1] = 37
pop_2[2,2] = 367
#Total rows and columns
pop_2[1,3] = pop_2[1,1] + pop_2[1,2]
pop_2[2,3] = pop_2[2,1] + pop_2[2,2]
pop_2[3,1] = pop_2[1,1] + pop_2[2,1]
pop_2[3,2] = pop_2[1,2] + pop_2[2,2]
N_2 = sum(pop_2[1,1] + pop_2[1,2] + pop_2[2,1] + pop_2[2,2])
pop_2[3,3] = N_2
pop_2
##                  Tine_Test_Pos Tine_Test_Neg Total
## Mantoux_Test_Pos           887            31   918
## Mantoux_Test_Neg            37           367   404
## Total                      924           398  1322

Hui-Walter model

A particular model formulation that was originally designed for evaluating diagnostic tests in the absence of a gold standard
Also known as the two_test - two_population setting/paradigm

Model Specification (‘hw_definition’)

hw_definition <- c("model{
  Population_1 ~ dmulti(prob_1, N_1)
  Population_2 ~ dmulti(prob_2, N_2)
  
  #Population_1
  
  # Test1+ Test2+
    prob_1[1] <- (prev[1] * ((se[1])*(se[2]))) + ((1-prev[1]) * ((1-sp[1])*(1-sp[2])))
  
  # Test1+ Test2-
    prob_1[2] <- (prev[1] * ((se[1])*(1-se[2]))) + ((1-prev[1]) * ((1-sp[1])*(sp[2])))

  # Test1- Test2+
    prob_1[3] <- (prev[1] * ((1-se[1])*(se[2]))) + ((1-prev[1]) * ((sp[1])*(1-sp[2])))

  # Test1- Test2-
    prob_1[4] <- (prev[1] * ((1-se[1])*(1-se[2]))) + ((1-prev[1]) * ((sp[1])*(sp[2])))
    
    #Population_2
  
  # Test1+ Test2+
    prob_2[1] <- (prev[2] * ((se[1])*(se[2]))) + ((1-prev[2]) * ((1-sp[1])*(1-sp[2])))
  
  # Test1+ Test2-
    prob_2[2] <- (prev[2] * ((se[1])*(1-se[2]))) + ((1-prev[2]) * ((1-sp[1])*(sp[2])))

  # Test1- Test2+
    prob_2[3] <- (prev[2] * ((1-se[1])*(se[2]))) + ((1-prev[2]) * ((sp[1])*(1-sp[2])))

  # Test1- Test2-
    prob_2[4] <- (prev[2] * ((1-se[1])*(1-se[2]))) + ((1-prev[2]) * ((sp[1])*(sp[2])))

  prev[1] ~ dbeta(1, 1)
  prev[2] ~ dbeta(1, 1)
  
  se[1] ~ dbeta(1, 1)I(1-sp[1], )
  sp[1] ~ dbeta(1, 1)
  se[2] ~ dbeta(1, 1)I(1-sp[2], )
  sp[2] ~ dbeta(1, 1)

  #data# Population_1, Population_2, N_1, N_2
  #monitor# prev, prob_1, prob_2, se, sp
  #inits# prev, se, sp
  }
  ")

Let’s run the model

library('runjags')

Population_1 <- as.numeric(pop_1[1:2,1:2])
Population_2 <- as.numeric(pop_2[1:2,1:2])


prev <- list(chain1=c(0.05,0.99), chain2=c(0.95,0.05))
se <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))
sp <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))

results <- run.jags(hw_definition, n.chains=2)

Remember to check convergence and effective sample size!

pt <- plot(results)

pt$`prev[1].plot1`

pt$`prev[1].plot3`

print(pt[["prev[1].plot1"]])

print(pt[["se[1].plot1"]])

print(pt[["sp[1].plot1"]])

print(pt[["sp[1].plot3"]])

summary(results)
##               Lower95      Median    Upper95        Mean          SD Mode
## prev[1]   0.015134996 0.028044913 0.04307840 0.028627279 0.007283481   NA
## prev[2]   0.690915330 0.715582840 0.74133190 0.715430590 0.012890610   NA
## prob_1[1] 0.014157200 0.026360418 0.04027954 0.026920260 0.006791926   NA
## prob_1[2] 0.007735796 0.017439528 0.02914599 0.018038765 0.005652170   NA
## prob_1[3] 0.002429108 0.008405407 0.01684339 0.008968083 0.003920497   NA
## prob_1[4] 0.927018653 0.946654656 0.96376544 0.946072892 0.009512237   NA
## prob_2[1] 0.643579183 0.669379993 0.69426946 0.669237567 0.012986677   NA
## prob_2[2] 0.019577493 0.028400701 0.03778528 0.028660257 0.004672813   NA
## prob_2[3] 0.015911645 0.023935304 0.03248369 0.024196034 0.004240767   NA
## prob_2[4] 0.253458655 0.277790209 0.30180763 0.277906142 0.012333790   NA
## se[1]     0.955584310 0.968754974 0.98063872 0.968429806 0.006411650   NA
## se[2]     0.951673572 0.966191131 0.98003042 0.965888714 0.007232896   NA
## sp[1]     0.970809468 0.982881713 0.99337004 0.982256434 0.005949219   NA
## sp[2]     0.983163647 0.992115396 0.99852033 0.991520346 0.004173473   NA
##                  MCerr MC%ofSD SSeff         AC.10      psrf
## prev[1]   6.671775e-05     0.9 11918 -0.0141690806 1.0001733
## prev[2]   1.216458e-04     0.9 11229 -0.0055674909 1.0000695
## prob_1[1] 6.193725e-05     0.9 12025 -0.0148276346 1.0001792
## prob_1[2] 6.160316e-05     1.1  8418 -0.0028020366 0.9999652
## prob_1[3] 4.231125e-05     1.1  8586  0.0024832450 1.0000058
## prob_1[4] 9.451507e-05     1.0 10129 -0.0185115477 1.0000470
## prob_2[1] 1.144714e-04     0.9 12871  0.0008074865 1.0000947
## prob_2[2] 4.372839e-05     0.9 11419 -0.0020506095 1.0000674
## prob_2[3] 3.980346e-05     0.9 11351  0.0143988167 1.0000089
## prob_2[4] 1.124378e-04     0.9 12033 -0.0017992025 1.0001151
## se[1]     6.565136e-05     1.0  9538  0.0190230422 0.9999648
## se[2]     7.659965e-05     1.1  8916 -0.0020043666 1.0000175
## sp[1]     6.558031e-05     1.1  8229 -0.0026330896 0.9999660
## sp[2]     4.569984e-05     1.1  8340  0.0043092761 0.9999978

Exercise 1

Run the hw_definition model under the following different scenarios and interpret the results in each case.
Change the priors for Se[1] and Sp[1] and try Beta(5,1).

Solution - Model Specification (‘hw_definition’)

hw_definition_1 <- c("model{
  Population_1 ~ dmulti(prob_1, N_1)
  Population_2 ~ dmulti(prob_2, N_2)
  
  #Population_1
  
  # Test1+ Test2+
    prob_1[1] <- (prev[1] * ((se[1])*(se[2]))) + ((1-prev[1]) * ((1-sp[1])*(1-sp[2])))
  
  # Test1+ Test2-
    prob_1[2] <- (prev[1] * ((se[1])*(1-se[2]))) + ((1-prev[1]) * ((1-sp[1])*(sp[2])))

  # Test1- Test2+
    prob_1[3] <- (prev[1] * ((1-se[1])*(se[2]))) + ((1-prev[1]) * ((sp[1])*(1-sp[2])))

  # Test1- Test2-
    prob_1[4] <- (prev[1] * ((1-se[1])*(1-se[2]))) + ((1-prev[1]) * ((sp[1])*(sp[2])))
    
    #Population_2
  
  # Test1+ Test2+
    prob_2[1] <- (prev[2] * ((se[1])*(se[2]))) + ((1-prev[2]) * ((1-sp[1])*(1-sp[2])))
  
  # Test1+ Test2-
    prob_2[2] <- (prev[2] * ((se[1])*(1-se[2]))) + ((1-prev[2]) * ((1-sp[1])*(sp[2])))

  # Test1- Test2+
    prob_2[3] <- (prev[2] * ((1-se[1])*(se[2]))) + ((1-prev[2]) * ((sp[1])*(1-sp[2])))

  # Test1- Test2-
    prob_2[4] <- (prev[2] * ((1-se[1])*(1-se[2]))) + ((1-prev[2]) * ((sp[1])*(sp[2])))

  prev[1] ~ dbeta(1, 1)
  prev[2] ~ dbeta(1, 1)
  
  se[1] ~ dbeta(5, 1)I(1-sp[1], )
  sp[1] ~ dbeta(5, 1)
  se[2] ~ dbeta(1, 1)I(1-sp[2], )
  sp[2] ~ dbeta(1, 1)

  #data# Population_1, Population_2, N_1, N_2
  #monitor# prev, prob_1, prob_2, se, sp
  #inits# prev, se, sp
  }
  ")

Population_1 <- as.numeric(pop_1[1:2,1:2])
Population_2 <- as.numeric(pop_2[1:2,1:2])


prev <- list(chain1=c(0.05,0.99), chain2=c(0.95,0.05))
se <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))
sp <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))

results_1 <- run.jags(hw_definition_1, n.chains=2)
# Remember to check convergence and effective sample size!
plot(results_1)  # Click backwards to view all plots

pt_1 <- plot(results_1)

pt_1$`prev[1].plot1`

pt_1$`prev[1].plot3`

print(pt_1[["prev[1].plot1"]])

print(pt_1[["se[1].plot1"]])

print(pt_1[["sp[1].plot1"]])

print(pt_1[["sp[1].plot3"]])

summary(results_1)
##               Lower95      Median    Upper95        Mean          SD Mode
## prev[1]   0.014826923 0.028134541 0.04306610 0.028749402 0.007438955   NA
## prev[2]   0.689752755 0.715504409 0.74021997 0.715359756 0.012947617   NA
## prob_1[1] 0.014773174 0.026470000 0.04116364 0.027035759 0.006934728   NA
## prob_1[2] 0.007476164 0.017319275 0.02904750 0.017931072 0.005663263   NA
## prob_1[3] 0.002048212 0.008427675 0.01675232 0.009013356 0.004000690   NA
## prob_1[4] 0.927204204 0.946435070 0.96494019 0.946019814 0.009646584   NA
## prob_2[1] 0.643645747 0.669390102 0.69435104 0.669222470 0.012932065   NA
## prob_2[2] 0.020076707 0.028527605 0.03795316 0.028734909 0.004622913   NA
## prob_2[3] 0.016287730 0.023785364 0.03256279 0.024051480 0.004207159   NA
## prob_2[4] 0.253883200 0.277833017 0.30226045 0.277991141 0.012376384   NA
## se[1]     0.955465800 0.969014218 0.98037901 0.968654657 0.006363729   NA
## se[2]     0.951839701 0.965934467 0.97960843 0.965740429 0.007153362   NA
## sp[1]     0.970736268 0.983012684 0.99345855 0.982375180 0.005961492   NA
## sp[2]     0.983262065 0.992083439 0.99895634 0.991469690 0.004255582   NA
##                  MCerr MC%ofSD SSeff         AC.10      psrf
## prev[1]   7.102587e-05     1.0 10970  0.0010270420 1.0000487
## prev[2]   1.220861e-04     0.9 11247  0.0081960861 1.0000494
## prob_1[1] 6.592711e-05     1.0 11064  0.0003184868 1.0000450
## prob_1[2] 6.215152e-05     1.1  8303  0.0154585463 0.9999592
## prob_1[3] 4.446067e-05     1.1  8097 -0.0029011841 1.0000610
## prob_1[4] 9.457490e-05     1.0 10404 -0.0097024767 1.0002971
## prob_2[1] 1.124436e-04     0.9 13227 -0.0012019476 1.0004216
## prob_2[2] 4.219356e-05     0.9 12004  0.0099645710 1.0003116
## prob_2[3] 4.173170e-05     1.0 10164 -0.0057853883 1.0004645
## prob_2[4] 1.117283e-04     0.9 12270  0.0046014984 1.0000558
## se[1]     6.870037e-05     1.1  8580 -0.0066740779 1.0005514
## se[2]     7.433448e-05     1.0  9261  0.0094238536 1.0002994
## sp[1]     6.615695e-05     1.1  8120  0.0153886562 0.9999518
## sp[2]     4.798811e-05     1.1  7864 -0.0019992100 1.0001384

Exercise 2

Remove the I(1-sp[1], ) and I(1-sp[2],) from the model and run it again. What happens now?

Solution Model Specification (‘hw_definition’)

hw_definition_2 <- c("model{
  Population_1 ~ dmulti(prob_1, N_1)
  Population_2 ~ dmulti(prob_2, N_2)
  
  #Population_1
  
  # Test1+ Test2+
    prob_1[1] <- (prev[1] * ((se[1])*(se[2]))) + ((1-prev[1]) * ((1-sp[1])*(1-sp[2])))
  
  # Test1+ Test2-
    prob_1[2] <- (prev[1] * ((se[1])*(1-se[2]))) + ((1-prev[1]) * ((1-sp[1])*(sp[2])))

  # Test1- Test2+
    prob_1[3] <- (prev[1] * ((1-se[1])*(se[2]))) + ((1-prev[1]) * ((sp[1])*(1-sp[2])))

  # Test1- Test2-
    prob_1[4] <- (prev[1] * ((1-se[1])*(1-se[2]))) + ((1-prev[1]) * ((sp[1])*(sp[2])))
    
    #Population_2
  
  # Test1+ Test2+
    prob_2[1] <- (prev[2] * ((se[1])*(se[2]))) + ((1-prev[2]) * ((1-sp[1])*(1-sp[2])))
  
  # Test1+ Test2-
    prob_2[2] <- (prev[2] * ((se[1])*(1-se[2]))) + ((1-prev[2]) * ((1-sp[1])*(sp[2])))

  # Test1- Test2+
    prob_2[3] <- (prev[2] * ((1-se[1])*(se[2]))) + ((1-prev[2]) * ((sp[1])*(1-sp[2])))

  # Test1- Test2-
    prob_2[4] <- (prev[2] * ((1-se[1])*(1-se[2]))) + ((1-prev[2]) * ((sp[1])*(sp[2])))

  prev[1] ~ dbeta(1, 1)
  prev[2] ~ dbeta(1, 1)
  
  se[1] ~ dbeta(1, 1)
  sp[1] ~ dbeta(1, 1)
  se[2] ~ dbeta(1, 1)
  sp[2] ~ dbeta(1, 1)

  #data# Population_1, Population_2, N_1, N_2
  #monitor# prev, prob_1, prob_2, se, sp
  #inits# prev, se, sp
  }
  ")

Population_1 <- as.numeric(pop_1[1:2,1:2])
Population_2 <- as.numeric(pop_2[1:2,1:2])


prev <- list(chain1=c(0.05,0.99), chain2=c(0.95,0.05))
se <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))
sp <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))

results_2 <- run.jags(hw_definition_2, n.chains=2)
# Remember to check convergence and effective sample size!
plot(results_2)  # Click backwards to view all plots

pt_2 <- plot(results_2)

pt_2$`prev[1].plot1`

pt_2$`prev[1].plot3`

print(pt_2[["prev[1].plot1"]])

print(pt_2[["se[1].plot1"]])

print(pt_2[["sp[1].plot1"]])

print(pt_2[["sp[1].plot3"]])

print(pt_2[["se[2].plot1"]])

summary(results_2)
##               Lower95     Median    Upper95        Mean          SD Mode
## prev[1]   0.015011314 0.02810646 0.04296564 0.028699435 0.007325016   NA
## prev[2]   0.690992901 0.71563005 0.74209896 0.715575887 0.013008858   NA
## prob_1[1] 0.014169211 0.02643493 0.04025525 0.026983627 0.006826925   NA
## prob_1[2] 0.007875020 0.01739199 0.02930559 0.017992973 0.005620845   NA
## prob_1[3] 0.002184289 0.00832246 0.01673989 0.008913238 0.003969325   NA
## prob_1[4] 0.927010407 0.94661370 0.96400615 0.946110162 0.009511643   NA
## prob_2[1] 0.643445447 0.66943086 0.69458122 0.669312627 0.013060432   NA
## prob_2[2] 0.019986687 0.02843649 0.03768666 0.028651309 0.004561164   NA
## prob_2[3] 0.016203085 0.02396447 0.03282358 0.024238713 0.004293344   NA
## prob_2[4] 0.253234837 0.27767702 0.30196209 0.277797351 0.012439240   NA
## se[1]     0.955313036 0.96868036 0.98048564 0.968349287 0.006465680   NA
## se[2]     0.952246961 0.96611973 0.97977093 0.965881641 0.007055120   NA
## sp[1]     0.970243662 0.98293571 0.99280968 0.982307275 0.005915084   NA
## sp[2]     0.983237189 0.99219836 0.99874894 0.991581913 0.004223711   NA
##                  MCerr MC%ofSD SSeff         AC.10      psrf
## prev[1]   7.051710e-05     1.0 10790  7.674084e-05 1.0003675
## prev[2]   1.249699e-04     1.0 10836 -3.169610e-04 1.0003569
## prob_1[1] 6.537476e-05     1.0 10905 -4.468137e-04 1.0003389
## prob_1[2] 6.259661e-05     1.1  8063 -2.809511e-04 1.0001406
## prob_1[3] 4.322508e-05     1.1  8433  7.642535e-03 1.0002174
## prob_1[4] 9.264131e-05     1.0 10541  5.740739e-04 1.0000597
## prob_2[1] 1.187534e-04     0.9 12095  4.492070e-03 1.0003792
## prob_2[2] 4.101602e-05     0.9 12366  1.028519e-03 0.9999808
## prob_2[3] 4.206362e-05     1.0 10418  2.047872e-02 1.0000397
## prob_2[4] 1.149664e-04     0.9 11707  1.745897e-03 1.0002765
## se[1]     6.921380e-05     1.1  8727  2.335653e-02 1.0003048
## se[2]     7.266383e-05     1.0  9427  7.800848e-03 0.9999859
## sp[1]     6.640189e-05     1.1  7935  1.804983e-04 1.0001219
## sp[2]     4.661818e-05     1.1  8209  7.924786e-03 1.0003453

Exercise 3

Try to run the model with different initial values, that explore the whole parameter space and removing I(1-sp[1],).

For example try it with:

se <- list(chain1=c(0.01,0.99), chain2=c(0.99,0.01))
sp <- list(chain1=c(0.01,0.99), chain2=c(0.99,0.01))

Solution

The model is the same like as in the above scenario (hw_definition_2)

Population_1 <- as.numeric(pop_1[1:2,1:2])
Population_2 <- as.numeric(pop_2[1:2,1:2])


prev <- list(chain1=c(0.05,0.99), chain2=c(0.95,0.05))
se <- list(chain1=c(0.01,0.99), chain2=c(0.99,0.01))
sp <- list(chain1=c(0.01,0.99), chain2=c(0.99,0.01))

results_3 <- run.jags(hw_definition_2, n.chains=2, sample = 100000)
# Remember to check convergence and effective sample size!
plot(results_3)  # Click backwards to view all plots

pt_3 <- plot(results_3)

pt_3$`prev[1].plot1`

pt_3$`prev[1].plot3`

print(pt_3[["prev[1].plot1"]])

print(pt_3[["se[1].plot1"]])

print(pt_3[["sp[1].plot1"]])

print(pt_3[["sp[1].plot3"]])

print(pt_3[["se[2].plot1"]])

summary(results_3)
##               Lower95      Median    Upper95        Mean          SD Mode
## prev[1]   0.015154582 0.028079649 0.04294688 0.028632057 0.007270327   NA
## prev[2]   0.690483159 0.715630958 0.74086809 0.715340989 0.012949846   NA
## prob_1[1] 0.014662730 0.026397736 0.04058197 0.026919481 0.006779520   NA
## prob_1[2] 0.007950100 0.017411169 0.02927751 0.018002716 0.005609978   NA
## prob_1[3] 0.002160442 0.008386046 0.01666301 0.008947818 0.003947023   NA
## prob_1[4] 0.927542366 0.946665368 0.96424735 0.946129985 0.009514629   NA
## prob_2[1] 0.642836937 0.669039765 0.69347360 0.669034858 0.012993111   NA
## prob_2[2] 0.020017827 0.028496244 0.03771279 0.028727058 0.004575018   NA
## prob_2[3] 0.015947471 0.023963876 0.03237984 0.024224011 0.004216653   NA
## prob_2[4] 0.253945453 0.277856590 0.30219512 0.278014073 0.012380665   NA
## se[1]     0.956013612 0.968677497 0.98071913 0.968372692 0.006350692   NA
## se[2]     0.951911789 0.965938232 0.97953781 0.965773842 0.007061216   NA
## sp[1]     0.970692441 0.982929910 0.99304449 0.982298017 0.005902372   NA
## sp[2]     0.983253015 0.992136762 0.99871786 0.991542641 0.004195551   NA
##                  MCerr MC%ofSD SSeff         AC.10      psrf
## prev[1]   5.140897e-05     0.7 20000  2.851607e-04 1.0004309
## prev[2]   9.222759e-05     0.7 19715  3.357144e-03 1.0000006
## prob_1[1] 4.793845e-05     0.7 20000  2.100863e-04 1.0004153
## prob_1[2] 4.006850e-05     0.7 19603  1.397071e-03 0.9999973
## prob_1[3] 2.814228e-05     0.7 19671 -1.038447e-03 1.0001594
## prob_1[4] 6.727859e-05     0.7 20000  2.223914e-04 1.0002099
## prob_2[1] 9.187517e-05     0.7 20000 -5.566153e-04 1.0000736
## prob_2[2] 3.235026e-05     0.7 20000  8.991125e-04 0.9999931
## prob_2[3] 2.981624e-05     0.7 20000 -8.326157e-04 1.0003006
## prob_2[4] 8.754452e-05     0.7 20000  2.277412e-03 1.0000170
## se[1]     4.490618e-05     0.7 20000 -3.803709e-05 1.0002579
## se[2]     4.993034e-05     0.7 20000  4.214278e-03 0.9999810
## sp[1]     4.214329e-05     0.7 19615  1.900281e-03 1.0000090
## sp[2]     2.992987e-05     0.7 19650 -8.217058e-04 1.0003295

The model fails to converge to one solution, but reports two solution that are complementary

Exercise 4

Run the model with only 1 population (either pop_1 or pop_2). What happens then?

We’ll remove population 2

Solution

Model Specification (‘hw_definition’)

hw_definition_4 <- c("model{
  Population_1 ~ dmulti(prob_1, N_1)

  #Population_1
  
  # Test1+ Test2+
    prob_1[1] <- (prev[1] * ((se[1])*(se[2]))) + ((1-prev[1]) * ((1-sp[1])*(1-sp[2])))
  
  # Test1+ Test2-
    prob_1[2] <- (prev[1] * ((se[1])*(1-se[2]))) + ((1-prev[1]) * ((1-sp[1])*(sp[2])))

  # Test1- Test2+
    prob_1[3] <- (prev[1] * ((1-se[1])*(se[2]))) + ((1-prev[1]) * ((sp[1])*(1-sp[2])))

  # Test1- Test2-
    prob_1[4] <- (prev[1] * ((1-se[1])*(1-se[2]))) + ((1-prev[1]) * ((sp[1])*(sp[2])))
    

  prev[1] ~ dbeta(1, 1)
  se[1] ~ dbeta(1, 1)I(1-sp[1],)
  sp[1] ~ dbeta(1, 1)
  se[2] ~ dbeta(1, 1)I(1-sp[2],)
  sp[2] ~ dbeta(1, 1)

  #data# Population_1, N_1
  #monitor# prev, prob_1,  se, sp
  #inits# prev, se, sp
  }
  ")

Population_1 <- as.numeric(pop_1[1:2,1:2])


prev <- list(chain1=c(0.05), chain2=c(0.95))
se <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))
sp <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))

results_4 <- run.jags(hw_definition_4, n.chains=2)
# Remember to check convergence and effective sample size!
plot(results_4)  # Click backwards to view all plots

pt_4 <- plot(results_4)

pt_4$`prev.plot1`

pt_4$`prev.plot3`

print(pt_4[["prev.plot1"]])

print(pt_4[["se[1].plot1"]])

print(pt_4[["sp[1].plot1"]])

print(pt_4[["sp[1].plot3"]])

print(pt_4[["se[2].plot3"]])

print(pt_4[["se[2].plot4"]])

print(pt_4[["se[2].plot5"]])

summary(results_4)
##               Lower95      Median    Upper95       Mean          SD Mode
## prev      0.019013278 0.038341165 0.06161965 0.03946337 0.011228794   NA
## prob_1[1] 0.012043887 0.023475571 0.03662748 0.02406689 0.006416080   NA
## prob_1[2] 0.008585623 0.018714963 0.03110883 0.01933945 0.005895514   NA
## prob_1[3] 0.003211859 0.009860699 0.01886579 0.01044148 0.004226694   NA
## prob_1[4] 0.927745026 0.946748161 0.96475071 0.94615218 0.009599461   NA
## se[1]     0.634003775 0.852750996 0.99999306 0.83821366 0.109583867   NA
## se[2]     0.506944192 0.754465929 0.99999240 0.75081035 0.144351054   NA
## sp[1]     0.976337501 0.989594679 0.99999630 0.98892152 0.007029213   NA
## sp[2]     0.986310908 0.994833454 0.99999952 0.99411332 0.004161227   NA
##                  MCerr MC%ofSD SSeff        AC.10      psrf
## prev      1.661421e-04     1.5  4568  0.022368201 1.0003222
## prob_1[1] 5.537503e-05     0.9 13425 -0.013294259 1.0004676
## prob_1[2] 5.279877e-05     0.9 12468  0.002956750 1.0000889
## prob_1[3] 4.239210e-05     1.0  9941  0.002543146 1.0001484
## prob_1[4] 8.900707e-05     0.9 11632 -0.001391273 1.0000299
## se[1]     1.471751e-03     1.3  5544 -0.008164718 0.9999978
## se[2]     2.344950e-03     1.6  3789  0.045963785 1.0002313
## sp[1]     1.097099e-04     1.6  4105  0.041163880 1.0001345
## sp[2]     5.603503e-05     1.3  5515 -0.002938294 1.0000742

How do the results look? # Compare the results with the scenarios above? Try also by removing population 1